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Abstract  -  A  novel  segmentation  technique  which  may  be  useful  for 
two  dimensional  (2D)  magnetic  resonance  (MR)  image 
segmentation  is  presented.  The  technique  utilizes  a  dynamic  target 
tracking  algorithm  and  a  Kalman  filter  and  permits  edges  to  be 
followed  in  the  presence  of  intensity  variation  similar  to  that 
found  in  MR  images.  Segmentation  of  two  synthetic  test  images, 
one  with  intensity  nonuniformity  and  one  without,  is  performed. 
Fuzzy  c-means  clustering  with  pixel  intensity  features  is  used  to 
segment  the  same  test  images  for  qualitative  comparison. 

Keywords  -  Image  segmentation,  edge  tracing,  Kalman  filter, 
intensity  nonuniformity. 

I.  Introduction 

Magnetic  Resonance  (MR)  images  are  excellent  sources  of 
patient-specific  anatomical  information.  Automatic 
segmentation  of  these  images  into  component  tissue  classes 
provides  a  method  for  reproducible  extraction  of  this 
information.  One  problem  that  complicates  this  process, 
however,  is  intensity  nonuniformity,  an  artifact  in  MR  images 
which  is  evident  as  a  gradual  variation  in  intensity  over 
otherwise  identical  tissue  classes.  Intensity  nonuniformity  has 
several  causes,  notably,  inhomogeneity  in  radio  frequency  (RF) 
transmitter  and  receiver  coils  during  image  acquisition  [1], 

MR  images  provide  excellent  soft  tissue  contrast  so  that 
intensity-related  features  are  natural  choices  for  use  with 
automatic  segmentation  methods.  However,  compensation  for 
intensity  nonuniformity  must  be  included  in  order  for  such 
methods  to  be  effective. 

Although  it  is  possible  to  perform  some  compensation 
during  image  acquisition,  equipment  or  protocol  modifications 
are  typically  required.  Furthermore,  retrospective  application  of 
these  corrective  measures  is  not  possible.  Therefore, 
compensation  applied  as  a  post-processing  step  is  considered  to 
be  desirable  [2]. 

Adaptive  fuzzy  c-means  [2],  [3],  and  statistically-based 
methods  [4],  [5]  are  examples  of  techniques  which  have  been 
developed  to  perform  automatic  image  segmentation  in  the 
presence  of  MR  intensity  nonuniformity.  Other  methods,  such 
as  nonlinear  filtering  [6],  are  intended  to  address  the 
nonuniformity  independently,  permitting  subsequent 
segmentation  of  the  intensity  corrected  image. 

Image  segmentation  can  be  performed  by  voxel  labeling, 
involving  classification  of  each  image  voxel  or  by 
identification  of  the  bounding  surfaces  of  objects  in  the  image. 
The  adaptive  fuzzy  c-means  methods  [2],  [3]  and  the  statistical 
methods  [4],  [5]  are  examples  of  techniques  which  perform 
voxel  labeling. 

Determining  the  object  boundaries  in  two  dimensional 
images  can  be  done  by  application  of  active  contours  [10]  or  by 
edge  tracing  [11],  We  describe  a  technique  for  edge  tracing 


a)  Unbiased  image  b)  Gain  Field 


c)  Biased  image 


Fig.  1  Synthetic  test  images. 

which  includes  a  Kalman  filter  and  a  dynamic  target  tracking 
algorithm  to  associate  edge  pixels  into  object  boundaries. 

II.  Methodology 

A.  Synthetic  Images 

Fig.  1  panels  (a),  and  (c)  show  the  two  synthetic  test 
images.  The  shapes  of  the  objects  in  the  images  have  been 
chosen  to  resemble  cortical  gray  matter  and  white  matter  in  MR 
images  of  the  brain.  Each  image  has  size  of  200x200  pixels 
with  256  gray  levels.  The  unbiased  image  was  formed  by 
interpolating  a  small  set  of  points  with  cubic  splines  to  form 
boundaries  of  closed  regions.  These  boundaries  were  then 
converted  to  discrete  pixels  and  the  enclosed  regions  were 
filled  with  a  selected  gray  level  value.  The  unbiased  test  image 
has  three  gray  levels  with  a  difference  of  fifty  gray  levels 
between  the  brightest  region  and  the  intermediate  intensity 
level. 

In  MR  images,  intensity  nonuniformity  has  been 
approximated  by  exponential  functions  [1],  For  the  test  images, 
the  gain  field  ( g ),  which  simulates  the  intensity 
nonuniformity,  was  formed  using  a  two  dimensional 
exponential  function: 

g  =  ^V{-kl{X-Xc)2-k2{Y-Yc)2)  (1) 

where  ( X ,  Y)  are  pixel  coordinates  and  (Xc ,  Yc )  are  the 
coordinates  of  the  image  centre.  Parameters  /c  |  and  k2  were 
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chosen  to  provide  a  fifty  percent  intensity  reduction  at  the 
image  edges  along  the  principal  axes. 

The  biased  image  was  formed  by  multiplying  the  unbiased 
image  by  the  gain  field,  thus  simulating  an  image  with  intensity 
nonuniformity. 

B.  Fuzzy  c-Means  Clustering 

Fuzzy  c-means  clustering  is  a  pattern  recognition  technique 
which  is  used  in  image  segmentation  [7].  Each  pixel  is 
evaluated  according  to  a  selected  feature  set,  forming  a  set  of 
vectors,  or  a  set  of  points,  in  the  feature  space.  Clusters  are 
regions  in  the  feature  space  with  a  high  density  of  such  points. 
A  prototype  vector,  or  cluster  centre  for  each  cluster  is  found 
by  an  iterative  computation  that  minimizes  the  objective 
function 

F=iiu™xDik2  (2) 

k=U= 4 

where  c  is  the  number  of  clusters  to  form,  N  is  the  number  of 
pixels  in  the  image,  m  is  the  fuzzification  factor  (typically 
chosen  to  be  2),  uik  is  the  membership  of  pixel  k  in  cluster  i , 

Dik  =  |  fk  -  Vj  |  is  the  distance  between  the  k  th  feature  vector 

(  fk )  and  the  /  th  prototype  vector  (v, ) .  Once  the  prototype 
vectors  and  membership  values  have  been  found,  each  pixel 
can  be  assigned  to  the  class  of  maximum  membership  to 
complete  the  segmentation. 

Fuzzy  c-means  clustering  is  used  here  to  demonstrate  the 
problem  that  can  occur  when  image  segmentation  is  performed 
without  attention  to  the  effect  of  nonuniform  intensity.  The 
features  selected  for  each  pixel  are  the  pixel  intensity  and  the 
intensity  of  the  four  nearest  neighbours.  The  number  of 
clusters  is  three  (ie.c  =  3)  and  the  fuzzification  factor  is  2. 

C.  Dynamic  Edge  Tracing 

Segmentation  by  edge  tracing  involves  edge  detection 
followed  by  association  of  edge  pixels  into  object  boundaries. 
In  the  case  of  the  edge  tracing  method  described  here,  edge 
detection  is  performed  on  the  input  image  followed  by  a  line  by 


Fig.  2.  Tracking  System  Block  Diagram 


line  scan  of  the  edge  image.  On  each  line  scan,  the  edge  pixel 
positions  and  image  intensity  at  the  edge  pixel  are  used  as  input 
to  a  multiple  target,  dynamic  tracking  algorithm.  The  scanning 
procedure  introduces  a  history,  allowing  edges  in  the  image  to 
be  followed  along  what  amounts  to  a  time  dimension.  Each 
new  line  brings  a  set  of  updated  edge  positions  which  are  then 
associated  with  existing  edge  data  from  the  previous  lines. 
Edge  positions  that  cannot  be  associated  with  an  existing  track 
are  used  to  start  two  new  tracks,  one  to  the  left  and  one  to  the 
right  of  the  scan  direction.  Tracks  which  follow  the  same  edge 
but  in  different  directions  will  terminate  on  each  other.  These 
occurrences  and  the  common  start  points  are  used  to  assemble 
the  edges  at  the  end  of  the  scan. 

The  effective  “movement”  of  the  edge  from  one  line  to  the 
next  during  the  image  scan  simulates  a  dynamic  system. 
Dynamic  systems,  linear  or  nonlinear,  are  described  by  state 
and  state  transitions.  State  is  a  quantitative  description  of  past 
and  present  behaviour,  sufficient  information  to  predict  future 
behaviour.  State  transition  is  the  description  of  how  one  state  is 
transformed  into  another.  For  example,  in  aircraft  tracking, 
state  would  include  the  position  and  velocity  of  the  aircraft. 
The  aircraft  position  could  be  predicted  at  a  future  time  based 
on  its  current  position  and  current  velocity. 

Automatic  tracking  algorithms  are  normally  used  to  monitor 
the  movement  of  aircraft  or  other  targets  of  interest  [8].  In  the 
classical  target  tracking  application,  sensors  provide 
measurements  of  the  target  state  (eg.  position  and  velocity)  to 
the  tracking  system  at  equal  time  intervals.  The  measured  target 
data  is  compared  to  predicted  target  data  and  if  sufficient 
correlation  exists,  the  measured  data  is  incorporated  into  the 
target  history  and  a  new  prediction  is  formed  for  the  next  input 
sample.  In  this  way,  observations  taken  at  different  times  can 
be  associated  together  and  the  path  taken  by  the  target  can  be 
followed. 

The  functions  of  the  tracking  system  are  data  association 
and  state  estimation.  Data  association  is  the  process  by  which 
new  data  is  correlated  with  existing  data  and  the  path  of  a  target 
is  updated.  State  estimation  is  the  process  whereby  a  target 
state  estimate  is  computed  using  a  priori  noise  statistics  and 
past  samples  and  whereby  a  predicted  target  state  is 
determined.  The  target  state  estimate  and  next  sample 
prediction  are  produced  by  a  tracking  filter  with  the  prediction 
presented  to  the  data  association  process  at  the  next  time 
interval.  A  block  diagram  of  such  a  tracking  system  is  shown  in 
Fig.  2. 

In  the  edge  tracing  method  described  here,  the  tracking 
filter  used  for  state  estimation  is  a  Kalman  filter  [9],  the 
recursive  solution  to  the  discrete  time,  linear,  minimum 
variance  estimation  problem  and  the  statistical  estimator  most 
often  used  in  dynamic  tracking  [8].  For  a  given  track,  the 
Kalman  filter  is  used  to  predict  the  edge  position  on  the  next 
line,  facilitating  the  association  of  the  next  set  of  edge  positions 
into  the  existing  edge  tracks. 

The  Kalman  filter  is  defined  with  the  assumptions  of  a 
linear,  dynamic  system  and  zero  mean,  gaussian  noise. 
Gaussian  distributions  and  linear  dynamics  are  natural 
assumptions  especially  if  statistical  data  is  not  largely  available 
[9], 
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The  Kalman  filter  can  be  used  to  estimate  the  state  of  a 
discrete  process  that  is  governed  by  the  linear,  stochastic 
difference  equation 

Xk+l=Akxk+wk  (3) 

where  k  is  the  step  counter,  x k  is  the  state  vector  at  step  k , 
Ak  is  the  state  transition  matrix,  and  wk  is  the  process  noise 
vector.  The  process  noise  is  assumed  to  be  zero  mean  with 
gaussian  statistics.  Measurements  related  to  the  target  are  also 
assumed  to  contain  zero  mean,  gaussian  noise: 

zk=Hkxk+vk  (4) 

where  zk  is  the  measurement  vector  at  step  k ,  II  k  is  the 
measurement  matrix,  and  vk  is  the  measurement  noise  vector. 

The  actual  state  of  the  target  is  not  known  and  the  Kalman 
filter  is  used  to  estimate  it  from  the  measurement  and  a 
previously  determined  state  prediction.  The  estimate  is  taken  to 
be  a  linear  combination  of  the  prediction  and  the  difference 
between  the  measurement  and  the  predicted  measurement. 

xk  =  xk  +  Kk  (zk  -  H kxk  )  (5) 

where  xk  is  the  state  estimate  at  step  k ,  xk  is  the  state 
prediction  at  step  k ,  and  Kk  is  the  Kalman  filter  gain.  Also, 
the  error  covariance  matrices  are  given  by: 

pk  =E{(xk  ~xk)(xk  - xk)T }  (6) 

pk  =E{{xk-xk\xk-xk)T)  (7) 

where  Pk  is  the  a  posteriori  error  covariance  matrix,  Pk  is  the 
a  priori  error  covariance  matrix,  and  E{}  represents 
mathematical  expectation.  The  Kalman  filter  gain  ( Kk )  is 
determined  by  minimization  of  the  error  covariance  matrix 

(pk  )  [9],  [12]. 

At  each  measurement  interval  the  Kalman  filter  gain  matrix, 


state  estimate  vector,  and  error  covariance  matrix  are  updated. 

Kk  =  pk~HkT  {HkPkHkT  +  Rk  )-1  (8) 

xk=xk  +Kk(zk-Hkxk)  (9) 

Pk=(I-KkHk)Pk  (10) 

where  Rk  =E{vkvk  }  is  the  measurement  noise  covariance 
matrix  and  /  is  the  identity  matrix. 

Prediction  of  the  next  state  is  also  done  at  each  time  interval 

xk+\=Akxk  (n) 

Pk+ 1  =  AkPkAkT  +Qk  (12) 


T 

where  Qk=E{wkwk  }  is  the  process  noise  covariance 
matrix. 

The  data  association  process  will  utilize  the  predicted  error 
covariance  matrix  to  form  a  bounding  window  around  the 
predicted  measurement.  Any  measurement  that  appears  within 
this  window  is  a  candidate  for  association. 

We  use  a  simple  two  state  filter  where  target  state  consists 
of  position  and  velocity  and  the  measurement  is  of  position 


a)  Edge  tracing 

-  unbiased  image. 


c)  Edge  tracing 
-  biased  image. 


b)  Fuzzy  clustering 
-  unbiased  image. 


d)  Fuzzy  clustering 
-  biased  image. 


Fig.  3.  Segmentation  results. 

only.  Under  these  conditions,  Ak  and  H k  can  easily  be 
defined.  That  is,  for 

xk=[X,Y,X,Y?  (13) 

where  (X,  Y )  represents  the  target  position  in  two  dimensions 

and  (X,  Y)  represents  the  target  velocity  in  two  dimensions 
and  assuming  a  unit  time  step,  the  state  transition  and 
measurement  matrices  become: 

'10  10 
0  10  1 

Ak  ~ 

0  0  10 
0  0  0  1 

which  is  to  say  that  the  next  target  state  will  be  estimated  from 
the  current  position  and  current  velocity,  that  only  position  is 
measured,  and  that  no  particular  measurement  correction  is 
required.  Furthermore,  these  two  matrices  will  remain  constant 
for  all  k . 

Our  calculations  are  done  in  this  manner  with  the  exception 
that  three  dimensions  are  used,  these  being  the  two  coordinates 
of  the  edge  pixel  and  the  image  intensity  at  the  edge  pixel 
location. 


(14),  Hk  = 


10  0  0 
0  10  0 


(15) 


III.  Results 

Fig.  3  shows  the  results  from  segmentation  of  the  unbiased 
and  biased  test  images  by  the  edge  tracing  technique  and  by  the 
fizzy  clustering  approach.  Panels  (b)  and  (d)  show  that 
clustering  works  well  when  the  intensity  is  uniform  (b)  but  that 
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given  sufficient  intensity  nonuniformity,  errors  occur  in  the 
pixel  assignments.  In  Panel  (d),  peripheral  portions  of  the  high 
intensity  region  are  classed  with  lower  intensity  pixels  and  the 
central  portion  of  the  high  intensity  region  is  expanded, 
exhibiting  a  circularly  shaped  artifact  due  to  the  exponentially 
shaped  intensity  variation. 

The  results  from  the  dynamic  edge  tracing  algorithm  are 
shown  in  (a)  and  (c).  In  each  case,  the  high  intensity  region  is 
outlined  with  a  black  contour  and  the  medium  intensity  region 
is  outlined  with  a  white  contour.  These  lines  coincide  with  the 
edge  pixels  very  well. 

The  fit  was  evaluated  by  reconstructing  the  test  image  using 
each  of  the  two  sets  of  edge  contours.  Upon  comparison  with 
the  original,  unbiased  test  image,  it  was  found  that  the 
reconstruction  using  contours  from  the  unbiased  image 
segmentation  contained  one  pixel  classified  incorrectly.  The 
reconstruction  using  contours  from  the  biased  image 
segmentation  contained  ten  misclassified  pixels,  amounting  to 
0.025  percent  of  all  image  pixels. 

IV.  Discussion 

Quantitative  comparison  of  the  fuzzy  clustering  result  with 
the  edge  tracing  result  is  not  intended,  however,  having  the  two 
side  by  side  gives  an  opportunity  to  examine  the  advantages  of 
each.  Although  fuzzy  clustering  requires  that  the  number  of 
classes  be  known  a  priori,  it  can  be  extended  to  perform 
segmentation  of  three  dimensional  (3D)  images  relatively 
easily.  The  edge  tracing  method  does  not  require  the  number  of 
tissue  classes  to  be  known  a  priori  but  is  not  as  easily  extended 
to  3D. 

When  2D  images  are  considered,  operation  without  a  priori 
knowledge  of  the  number  of  tissue  classes  is  a  big  advantage 
especially  if  the  goal  is  automatic  analysis  of  images  where 
pathology  may  be  involved.  Fuzzy  clustering  may  not  position 
a  cluster  prototype  so  as  to  identify  a  relatively  small  region  of 
distinct  intensity  in  an  image  when  there  are  much  larger 
numbers  of  pixels  in  other  intensity  groups.  Consideration  of 
the  objective  function  in  (2)  will  confirm  this.  Since  the 
objective  function  is  based  on  minimizing  the  sum  of  all 
distances,  a  small  but  distinct  group  in  the  feature  space  may 
not  have  enough  accumulated  distance  to  attract  a  cluster 
centre.  The  edge  tracing  technique,  however,  would  find  all 
regions  where  there  is  an  identifiable  edge. 

Edge  tracing  methods  require  some  degree  of  edge 
continuity  to  be  successful  [11].  The  edge  tracing  technique 
described  here  does  not  require  adjacent  pixel  connectivity. 
Since  the  edge  position  is  permitted  to  have  a  “velocity”,  the 
current  velocity  for  that  edge  will  determine  where  the 
algorithm  searches  for  the  next  edge  pixel. 

V.  Conclusion 

Edge  detection  followed  by  a  line  by  line  scan  of  the  edge 
image  simulates  a  dynamic  system  where  the  edge  position 
“moves”  as  the  scan  proceeds.  The  application  of  a  dynamic 
tracking  algorithm  allows  the  edge  to  be  followed,  permitting 
edge  pixels  to  be  associated  into  object  boundaries.  Edge 


tracing  performed  in  this  manner  and  using  three  dimensions 
(the  two  edge  pixel  position  coordinates  and  the  image 
intensity  at  the  edge  pixel  position)  appears  to  be  a  viable 
method  for  2D  image  segmentation  in  the  presence  of  image 
intensity  nonuniformity. 
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